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I> ; Abstract 

A three dimensional attractive Bose-Einstein Condensate (BEC) is expected to collapse, when 
the number of the particles N in the ground state or the interaction strength Ao exceeds a critical 

w>: . 

value. We study systems of different particle numbers and interaction strength and find that 
even if the overall ground state is collapsed there is a plethora of fragmented excited states that 



are still in the metastable region. Utilizing the configuration interaction expansion we determine 
the spectrum of the ground ('yrast') and excited many-body states with definite total angular 
momentum quantum numbers ^ L ^ N and —L sj Ml ^ L, and we find and examine states 
that survive the collapse. This opens up the possibility of realizing a metastable system with 

<-> : 

overcritical numbers of bosons in a ground state with angular momentum L ^ 0. The multi-orbital 

mean-field theory predictions about the existence of fragmented metastable states with overcritical 

| numbers of bosons are verified and elucidated at the many-body level. The descriptions of the 

ON ' 

£f) • total angular momentum within the mean-field and the many-body approaches are compared. 
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I. INTRODUCTION 



Attractive trapped Bose-Einstein Condensates (BEC), since their first realization 




have gained increasing attention, due to the interesting phenomena they exhibit 
three dimensions a salient feature is the collapse of the condensate when the (negative) 
interacting energy per particle, is too large to be compensated by the kinetic energy The 
attraction brings the bosons so tightly close that the spatial extension of the wave function 



of the system shrinks to a point and the condensate eventually implodes 7H11I]. However, 
in trapped gases, metastable states, i.e., states that remain stable for some finite time, exist. 
The metastability of the condensate is a critical phenomenon: if at a given interaction 
strength A , the total number of the bosons exceeds a critical value N cr or vice versa (see 
for example 13]), the condensate will collapse. 

Much work has been devoted to exploring the ground metastable state of attractive 



BECs and its properties (see, for instance, [6|, Il4j-|l7j). Recent experiments have revealed 
new phenomena in attractive BECs that seem to go beyond the ground metastable state. 
In particular, it was found that states with over-critical number of bosons exist [lsj. It 
is natural to assume that excited states of the attractive BECs are involved. Furthermore, 
disagreements have been reported (see, e.g., fl) between the experiments and the predictions 
of the Gross-Pitaevskii (GP) theory on the critical value of the attraction strength, where the 
gas collapses. This motivates us to theoretically study excited states of attractive BECs. 
We go beyond the standard Mean-Field (MF) theory in the expectation of a thorougher 
understanding of the role of excited states in the stability of attractive BECs. 

To describe the statics and dynamics of a condensate one can adopt a mean-field (MF) 
approach, such as the famous Gross-Pitaevskii (GP) theory. The starting point of the GP 
description of a condensate is that all the bosons of the system reside in one and the same 
one-particle function (orbital) 0o( r ) an d hence the wave function $ of the system is merely 
a product of this prototypal orbital: 

$o(ri, r 2 , . . . r N ) = ^(r^n^) . . . </)q(t n ), (1) 

where N is the number of particles. The expectation value of the system's many-body (MB) 
Hamiltonian evaluated with this trial function is -E[$o] = (*&o\H\<&o)- The orbital that 
minizes this expectation value is the optimal orbital and is found to satisfy a non-linear, 



partial differential equation, the famous Gross-Pitaevskii equation 
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However, a GP ansatz is not, by definition, capable of describing fragmentation phenom- 



ena. The relaxation of the assumption that all bosons are condensed in 
found to lead to fruitful results: energetically favourable fragmented states 



r) has been 
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metastable states with overcritical number of bosons 



23j . excited 



24(, fermionized states and new Mott- 



insulator phases 25|, |26|]. Under this generalized MF approach the wave function of the 



system is rewritten as: 



$(ri, r 2 , . . . , r N ) = «S0i(ri)0i(r 2 ) . . . 0i(r ni )0 2 (r ni+ i) . . . 2 (r ni . 



Hi [FN), 



(2) 



where S is the symmetrizing operator. In such a description and in contrast to the GP 
case, M > 1 orbitals are allowed to be occupied by bosons. Within this multi-orbital Best 
Mean Field (BMF) approach J2I] the orbitals <pj of Eq. (j2j) as well as their occupations 
rij are calculated self- consistently, so as to minimize the total energy of the system. More 
specifically, the orbitals <pj are found to satisfy a system of M coupled non-linear differential 
equations, in place now of the GP equation. The wave function of the system [Eq. (j2J)], i.e., 
the symmetrized product of M different orbitals is sometimes called 'permanent', since it 



can be regarded as the permanent of a matrix (in direct analogy to the S 



27 



ater determinant). 
1 condensate, since 



A condensate, described by Eq. (J2J) is — generally — a fragmented 
the reduced one-particle density matrix will give 'signatures' in more than one eigenvalues. 
So, within the BMF framework fragmented states are well described and a GP state arises 
as an extreme case, where none but one orbital is occupied by all N bosons. Some relevant 
work to the BMF app roach, in favour of or against fragmentation of bosonic systems, can 



be found in Refs. 



23-122]. 



Going one step further, we write the wave function [ansatz) of the system as a linear 
combination \1/ = ^\ Cj$j of different states <J>j (permanents), which are taken from a 
set of orthonormal permanents $j, each one describing a condensate (fragmented or 
not) of N particles. This set which spans our configuration space, consists of all the 

permanents that result by distributing, in all possible ways, the N bosons over the M orbitals 
Such a state \1/ is known as a configuration interaction (CI) expansion 3^|. In contrast 
to the BMF, a CI state \l/ can further describe purely MB phenomena, such as depletion and 
fluctuations of the states. The coefficients Cj of the above expansion as well as the orbitals <j)j 
are determined variationally. Owing to the detailed analysis of Ref. 39j there is a formally 
defined Multi-Configurational Hartree method for Bosons (MCHB) and an efficient way to 
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determine the ground and excited states and their energies, i.e., the whole spectrum of a 
given Hamiltonian in a given configuration subspace of The expansion coefficients 



{Ci} are determined from the diagonalization of the respective secular matrix [38] while the 
one-particle wave functions 4>j are obtained by solving a system of M coupled non-linear 
(integro-) differential equations [39[]. The MCHB theory and its time-dependent counterpart 
have been successfully applied to a range of problems of one-dimensional ultracold boson 



gases, predicting various new phenomena 



4CH45| . However, it is, to this day, not feasible 



to exactly solve the MCHB equations in three-dimensional problems. To overcome this 
difficulty in the present work, as will be later explained, we implement a restricted version 
of the MCHB theory, namely a CI expansion of permanents, built over a one-parametric 
one-particle functions set. 

The complexity of a problem, in the framework of MCHB, depends firstly on the number 
M of the one-particle functions, that are to be determined variationally and secondly, on 
the symmetries of the system, which in general reduce the size of the configuration space. 
Since a complete configuration space would be infinite-dimensional, the space of the orbitals 
(and hence that of the permanents) has to be truncated and limited down to a relatively 
small number M, so that real calculations can be performed. Besides this limitation of 
the configuration space, we imply another constraint on the one-particle wave functions 
4>j\ we suppose that the one-particle states can be well approximated by wave functions 
— ansdtze — that are completely known, upon some real parameters Oi. In other words, 
we fix a priori the solutions of the system of M coupled non-linear differential equations 
to M mono-parametric families of complex wave functions and we then look for the values 
of <7j that 'optimize' the solution, i.e., values that assign an extremal value to the energy 



functional. These ansatze 
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46| are taken in our case after the (exact) solution of the 



corresponding non-interacting system, i.e., they are scaled Gaussians and inherit thus the 
symmetries of the original system. 

In an earlier relevant work 35[ Elgar0y and Pethick have derived and used a two-mode 
MB Hamiltonian, borrowed from the nuclear physics Lipkin model, to determine the ground 
state of an attractive trapped Bose gas. The modes correspond to the s-orbital loo an d 
the p-orbital Y" 10 and the Hamiltonian matrix is constructed over the set of permanents 
\n) = | no, N — no), where no bosons reside in the s-orbital and iV — no in the p-orbital, 
N being the total number of particles. Then, by rewriting the Hamiltonian in terms of a 
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quasispin operators J z , J + , j_ they calculated the population at each mode, in the ground 
state. All the configurations \n) are eigenstates of the quasispin operators with quantum 
numbers J = N/2 and J z = ~(iV — 2n ). The ground state, in the range of the parameters 
where it is not collapsed, was found to be not fragmented. However the authors did not 
examine excited states, which as we shall show in our work, can carry angular momentum 
(or 'non- minimum quasispin' in the case of 35 h and are metastable fragmented states that 
survive the collapse. Still, the work of Ref. \35\ has stimulated the present extended and 
more complete study. By including all the p-orbitals in our configuration space, we are 
able to write the wave function of the system as eigenfunction of the true total angular 
momentum operators and hence restore the symmetries of the problem. 

A set of relevant published works, where ultracold bosonic systems are examined with 
methods beyond the MF approach, includes Refs. 14H17j. However they pertain to (true 
or quasi-) two-dimensional systems, where the description of the angular momentum basis 
is fairly different and simpler than the analysis on a fully three-dimensional system that we 
present here. In addition, they do not examine the stability of the system with respect to 
the fragmentation of ground or excited states. 

To render our MB method more efficient we should take all the symmetries of the problem 
into account. The one-particle functions set M. that we use, i.e., the set of the cr^-orbitals, see 
Eqs. (l9])-f }TTj) below, consists of functions that have definite orbital angular momenta, mi and 
I, as well as parity (symmetry under spatial inversion). However one is more interested in the 
symmetries of the MB states since these will directly reduce the size of the configuration 
(Fock) space. 

We principally aim at investigating how the total angular momentum affects the stability 
and fragmentation of the system. To achieve this we first have to answer on what are the 
MB states with definite total angular momentum. We, therefore, define the MB operators 
L 2 ,L Z and their action on the permanents We then look for a {$j}-basis of MB states 
that are common eigenstates of L 2 ,L Z . Once the new basis is known we can rewrite 
the state \l/ and the Hamiltonian of the system on this basis for given eigenvalues L,Ml, 
i.e., over states with the same symmetry. In such a way the size of the new, rotated basis 
set {$f }, with the index L meaning hereafter that the members of this set have the same 
angular momentum quantum numbers, is significantly smaller than the original and 
the calculations are further facilitated (see also Appendix |A]). We will show that a general 
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state with definite angular momentum L,Ml, is a quantum MB state, i.e., a non MF 
state. 

We should also stress the relevance of the 'yrast' lines to the present work. The term 
'yrast' state (or level) has been coined to describe the lowest-in-energy states, for a given 
angular momentum, first in the context of nuclear physics and much later in the physics 
of ultracold Bose gases [48]. Herein we do explore the yrast states of attractive systems but 
we, also, look at the excited, i.e., non yrast states, for given L,Ml- However, the presence 
of attraction induces a subtle feature. Which state (ground, first excited, etc.) is accretided 
with the term 'yrast' will in principal change, as the absolute value of the interaction strength 
increases and the lowest-in-energy states start to collapse. 

An equally important goal is to examine the stability and the properties of the ground 
and excited states of different angular momenta of systems of various Ao and N. To do 
so, we employ the natural orbital analysis; the findings strongly support that states with 
angular momentum different than zero [large or not, depending on the quantity Xo(N — 1)] 
can exist, with a total number of particles well above the critical number of particles N^ p , 



as ca. 



culated from the GP theory. We verify, therefrom, the predictions of the BMF of 



Ref . 24J , that fragmented excited states exist and survive the collapse and we explain these 



features at the MB level. 

The structure of the paper is the following. In Sec. Ill Al we introduce our theoretical 
approach to (stationary) quantum bosonic gases and we define the MB and one-particle 
states that form the configuration spaces. In Sec. Ill Bl we give the expression for the total 
angular momentum operator, we derive a MB angular momentum basis, and we show how 
this partitions the configuration space. In Sec. Ill CI we define the main tools of the natural 
orbital analysis of the MB states. In Sec. II III the main results of this work, for systems 
of N = 12, 60 and iV = 120 are presented; in Sec. IIII Al MB states belonging to the same 
subspace L = are compared, while in Sec. IIII Bl we examine ground states of different 
L-subspaces. In Sec. IIII CI we further investigate the properties, namely fragmentation and 
variance of the expansion coefficients, see Eqs. (jSJ) and (|28|) below, of the previously found 
metastable MB states. In Sec. II VI we study the overall impact of the angular momentum on 
the state of the system with respect to its collapse, and we compare the role of the angular 
momentum within the MB and the MF theories. Last, Sec. IVl summarizes our results and 
provides concluding remarks. A set of relevant derivations are given in Appendices \M and 
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II. THEORY 



A. Preliminaries and basic definitions 



We consider a system of N identical spinless bosons of mass m confined by an exter- 
nal time-independent potential Vfa), and interacting with a general two-body interaction 
potential W(r$ — Tj), where i-j are the space coordinates of the i-th boson. 

The Hamiltonian of the system is: 

H = h + W, (3) 

with h = E?h(*i) = -£Ef V? 4 + Y?V{Ti) and W = TScjWin - r,) the many- 
body interaction operator. In the present work we choose V"(r$) = V(ri), i.e., the trap 
potential has spherical symmetry. For the interaction operator we will use the common 
delta function 5(r — r') representation, W(r — r') = Xo5(r — r'), where the parameter Ao 
measures the strength of the interparticle interaction. This parameter is proportional to the 
s-wave scattering length a s and takes on negative values for attractive interaction. Precisely, 
Ao = 47ra s y^|j, where oo is the frequency of the trapping potential. The time- independent 
MB Schrodinger equation reads: 

= EV, (4) 

where ^ = ^(ri, r 2 , . . . r^r) is the MB wave function of the system of N interacting bosons 
and E the eigenvalue of the operator H, corresponding to the state Even in the simple 
isotropic case analytic solutions for the MB Schrodinger equation are not known. Hence 
we will approximate the solution \1/ with a MB ansatz as an expansion over known states 
(permanents) 

N p 

\^) = j2am, (5) 

i 

where N p is the total number of permanents used in the expansion and Cj, % = 1 . . . N p the 
corresponding coefficients. The question that arises is what the set of MB basis functions 
consists of. Generally, this includes, as mentioned, all the permanents that result 
from distributing iV bosons over M orbitals. Readily, these states are those of Eq. (j2J). In 
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occupation number representation (Fock space representation) the same states take on the 
form: 

|$) = \n) = \n 1 ,n 2 , . . -,n M )- (6) 
Here rij denotes the respective occupation number of the one-particle functions (orbitals) 

Mathematically seen, the permanents |$) are the vectors that span the Fock space T of 
all iV-body wave functions. However it is possible to reduce the size of the 'working' Fock 
spaces, by partitioning the initial space into n- and L-subspaces, i.e., spaces of permanents 
with definite parity and angular momentum. The purpose to do so is twofold; firstly the 
resulting solution will possess the rotational symmetries of the system and secondly the 
working {$f } spaces are each time much smaller in size than the initial {$?} one. Note 
that, while this partioning of the configuration space with respect to angular momentum L 
is trivial in a two-dimensional system, it is not straightforward in the full 3D case, that we 
examine here. 

The MB Hamiltonian H in the {$j}-basis is represented as a secular matrix H, with 
elements: 

7Uj = (^\H\^). (7) 
By diagonalizing, i.e., solving the equation 

HQ = SC, (8) 

where C is the column vector of the expansion coefficient C = {C\, Ci . . . , Cn p } t and £ the 
respective eigenvalue, we obtain the energies (eigenvalues) and coefficients (eigenvectors) of 
the solutions of the system. Note that the Hamiltonian in Eq. ([7]) is evaluated still in the 
full basis of permanents. 

To complete the picture of our variational solutions, we give the one-particle function 
basis set, over which the permanents of Eq. (j2J) are constructed. This set of ansdtze consists 
of the known orbitals that solve the isotropic 3D quantum harmonic oscillator, scaled under 
a scaling parameter <ji. Precisely, it consists of four orbitals: the ground / = and the 
three I = 1 excited ones, which we scale with two parameters (a , ai), as have been already 



done in Ref. 24j. The parameters cr, will determine the shape (width) of the orbitals; 



their optimal values are such that, for a given set of coefficients C, the total energy of the 
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system takes an extremum. In this approximative way we restrict the solution of the system 
to functions that lie inside the monoparametric families of equations, which solve a scaled 
ordinary Schrodinger equation; the solution of a coupled system of nonlinear differential 
equations (MCHB equations) boils down to the determination of a set of parameters which 
minimizes the total energy. 

In this work much of the analysis of the quantities of the system (energy per particle, 
occupation numbers, variances of expansion coefficients) is done with respect to the scaling 
parameters. So for example, we expect to see a (local) minimum in the plot of E(ao, a{) 
against a , <j\ if the system is metastable, while the absence of minimum will signal a col- 
lapsing condensate. Moreover, the analysis of the resulting MB states (depletion, angular 
momentum) is performed at the optimal values of the parameters <7j, % = 0,1. 

The orbitals that we use, in Cartesian coordinates, have the form: 



(9) 



where 



and 



0i(r) = tp {x,ao)tpo(y,ao)(po(z,a ), 

H r ) = -^M^, <Tl)(Po(V, Vl)<Po(z, O-l ) + i<fo(x, (Tl)tfl(y, <Tl)(Po(z, ^l)), 

3 (r) = tp (x, (Ti)ip (y, cri)ipi(z, at), 
^(r) = -^=((^i(x,cri)v9o(l/,cx 1 )v9o(2,ai) - iip (x, <7i)</?i(y, (Ji)<fo(z, en)), 



Mx,°) = [—^) e (10) 



"4 


/mu\ 3 


7T 


\a 2 h) _ 



1/4 

xe 2 ^ (11) 



are orthonormal orbitals, i.e., (<fii\<fij) = i,j = 0,1. Here m is the particle mass and 
u is the trap frequency. Throughout this work the quantities used are dimensionless, i.e., 
h — m — u — 1. 

In terms of spherical harmonics Yi j7ni , i.e., under a change of coordinates, the orbitals of 
Eq. (ED are: 

<t> k {r)=<p l (r,<7 l )Y lmi (9,<t>), (12) 
where I = 0, 1, — I ^ mi ^ I and k = k(l, m{) = 1 + 1(1 + 1) — m/. 
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B. Angular momentum basis 



It is easy to see that the orbitals of Eq. (fl2l) constitute a set of common eigenstates 
of the orbital angular momentum operators L 2 ,L Z together with the parity (inversion) 
operator fl : II^(r) = r), with eigenvalues I = {0,1,1,1}, m ; = {0,1,0,-1} and 
7r = {1, —1, —1, —1}, respectively. 

We now want to express the total angular momentum operators at the MB level. For 
this purpose we switch to second quantization language and introduce the bosonic creation 
(annihilation) operators b\ (6j), associated with the orbital set {0j(r)} and which obey the 
usual bosonic commutation relations: b{bj — bjbi = 5ij. 

The total angular momentum operators are (see, e.g., j^): 

L 2 = L 2 Z + ^(L+L-+L-L + ), (13) 

Lz = J2 mibj lm b lrn n (14) 

l,mi 

L± = Y1 A{1 ' T™*)&L !± i 6 ^ , (15) 

l,mi 

where Ail, m\) = [(/ + mi)(l — m/ + l)] 1 ^ 2 and b\ m {bi mi ) creates (annihilates) a boson in the 
state 4>i mi , with orbital angular momentum quantum numbers I, mi. 
Applying Eq. f[T5|) to our basis of Eq. ([H]) with M = 4 we get: 

L 2 \n) = [n 2 (n 3 + 1) + ^3(^-4 + 1) + n 3 (n 2 + 1) + n 4 (n 3 + 1) + (n 2 - n 4 ) 2 ] |m, n 2 , n 3 , n 4 ) + 
+2^/713(713 - l)(n 2 + l)(n 4 + l)|ni, n 2 + 1, n 3 - 2, n 4 + 1} + 
J r1^n 2 n^(7%3 + l)(n 3 + 2~)\m,n 2 - l,n 3 + 2,n 4 - 1). 

(16) 

It can be easily seen that each permanent [Eq. ([6])] is an eigenstate of L z with eigenvalue 

M L = n 2 - n 4 , 

L x \n) = (n 2 -n 4 )|n). (17) 

But what happens to the eigenstates and eigenvalues of the L 2 operator? To answer this, 
one has to solve the eigenvalue equation: 

CC = AC, (18) 
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where C is the matrix representation of the operator L 2 in the basis of permanents [Eq. (J5j)] 
with matrix elements 

Cij = (fiil^lnj), (19) 

C is the column vector of the coefficients Ci and A the eigenvalue in question. 

A unitary transformation U will in general rotate the <3?-basis to a new <3> one. In this 
basis the secular matrix of Eq. ((7j) becomes: 

n itj = = j2 u lHk,iUi,i, (20) 

k,l 

where Uij are the matrix elements of U. If IA is simply the matrix of the eigenvectors of 
L then H, takes on the desired total-angular-momentum block-diagonal form. The above 
mentioned vector spaces, that the bases {$j}, {$j} span, are homomorphic and can be 
written as the direct sum of the subspaces {<E>f }: 

m = m = ©r l }. (2i) 

L 

We have numerically calculated the angular momentum states $j and the matrix transfor- 
mation U that transforms to the new basis $ = of eigenstates of L and Ml, for the cases 
of N = 12, 20, 60, 120 bosons. An analytic approach to the same problem of determining the 
states $j is presented in Appendix IB 11 For selected values of L and Ml we construct and 
diagonalize the block 1-L L of the secular Hamiltonian matrix % and find its eigenf unctions 
We use hereafter the index L to stress the fact that the state |\1/ L ) is an eigenstate of 
L 2 with quantum number L. The size of the block T-i L is found to be N p = N ~^ +2 (see also 
Appendix |A]), with the same number of eigenstates. We index the states with i, to 
denote the ground (i = 1) and the excited (1 < % < N ~^ +2 ) states belonging to this block of 
angular momentum L. When it is not transparent from the context, we will also use Aq cr or 
Ao :Cr , to denote the critical value of the interaction strength where the state ) collapses. 

C. Natural orbital analysis 

The first order reduced density matrix (RDM) for the state = ^n^™l^)> Eq. ([5]), is 
defined as: 

„ M 

p(r\r') = N / ^*(r\r 2 ...r N )^(r,r 2 ...r N )dr 2 dr 3 ...dr N = Y,Pv<l>i^')Hr) (22) 
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and the second order RDM: 

p(r 1; r 2 |ri, r' 2 ) = N(N -I) J **(ri, r' 2 , r 3 . . . r 7V )^(r 1 , r 2 , r 3 . . . r^drg ...dr N = 

M 

p lJ H0I(r , 1 )0*(r 2 )0 fc (r 1 )0 i (r 2 ), (23) 

with = (^\b\bj\^), pijki = (^\blbjb k bi\^J) being the elements of these matrices. The nat- 
ural orbitals 4>f° are defined as the eigenfunctions of p(r|r'), i.e., the one-particle functions 
that diagonalize the right hand side of Eq. (J22J) and their eigenvalues are known as natural 
occupations pi. In our system, the spherical symmetry of the Hamiltonian induces a zero 
"H-matrix element between states of different symmetry (angular momentum and parity). 
Hence p^ is diagonal and the ansatz orbital-set of Eq. (Q coincides with the set of the 
natural orbitals {(fif }, i = 1, • • • , M. 
The diagonal elements 

Pa = Pi = ^ CtCfjUj = (hi), (24) 

i — 1, . . . , M, are the expectation values of the number operators associated with the orbitals 
4>i and are, as explained, the natural occupations pi of the respective natural orbitals. We 
define the depletion d,i of the z-orbital as: 

di = l- pi/N. (25) 

The depletion di is an informative quantity which measures the relative number of particles 
that are depleted from the z-orbital. Throughout this work we extensively use the quantity 
di and also refer to it as the s-depletion. The diagonal elements 

Pun = J2 C n C M - rii) = (n, 2 ) - (fa) (26) 

n 

are related to the variance of the distribution of the coefficients C, of a given state |\&), 
associated with the orbitals 4>u as: 

T i = (^ 2 ) - (^i) 2 = Pan + Pui 1 ~ Pa)- ( 27 ) 

The norm 

M = sJtI + tI + tI + tI (28) 
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gives a measure of the variances [57| of all four orbitals of a state The r^'s and \r\ are 
highly useful quantities, as they measure the fluctuations around the occupation numbers of 



the natural orbita 
|r| is simply zero 



(t>f°. In the case of a mean- field state, where there are no fluctuations, 
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III. MANY-BODY RESULTS 

In this section we implement the many-body method described above, for systems of 
trapped ultra-cold gases. We present and discuss calculations regarding systems of iV = 
12, 60 and 120 bosons, embedded in a spherically symmetric trap. First, the interaction 
strength Aq is chosen each time such that the product |Ao|iV is kept fixed to the value 10.104 



59| | . This choice will permit a direct comparison of our results to those of Ref. [24] . where 
an attractive system of |Ao|./V = 10.104 was also examined. Later on, also other values of 
|Ao|iV are considered for the shake of completeness. In the following we examine states of 
definite angular momentum L, and positive parity II only. The latter makes the total 
angular momentum of each state increase at an even step, i.e., L = 0, 2, 4, . . . , N. 

A. Ground and excited states of the 'block' L=0 

For each of the above systems we examine states with definite angular momentum. We 
first calculate the energy per particle e = S/N of those states, as a function of the variational 
parameters o"o, (J\ of the orbitals [see Eqs. fflUj) and pip ]. We then look for the minimum eo 
of the energy with respect to these parameters. As mentioned earlier, the total absence of a 
minimum indicates unbound (total) energy and a collapsing system. Namely, as o"o, <J\ — > 
the orbitals of Eq. (Q contract to pointlike distributions. When existent, the minima are 
expected to be local only; an energy barrier separates the metastability from the collapse 
regions. The shape of the energy barrier determines the tunneling time of the system 
through this barrier and — generally — the higher the barrier is, the longer the system is 
expected to survive in this state. The variation of <Jq, <Ji takes place over states of the same 
symmetry and hence the surfaces ought not to cross (see Fig. [TJ See also the theoretical 



discussion on non-crossing of energy surfaces in 50|, [51] and references therein). Notice that, 



owing to the attractive interparticle interaction, the wave function of the system has to be 
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spatially shrunk, compared to that of the non-interacting system; indeed the optimal scaling 
parameters of the orbitals, are always found to obey ctq, <J\ < 1. 

The first system studied is that of N = 120 bosons, with attractive interaction of strength 
Ao=-0.0842. The energies per particle e(cr , ai) of three distinct states of this system are 
collectively presented in Fig. [TJ We first pick the state with quantum numbers L=0 and 
Ml=0 of the operators L 2 ,L Z , respectively. We find that the ground state is collapsed 
(lowest surface in Fig. [p. As the introduction of Sec. fl] suggests, we expect to find excited, 
fragmented states that can survive this collapse. Indeed, an examination of the spectrum 
of the states of the Hamiltonian of Eq. (j7|) reveals that the l^^o) excited state is the first 
to demonstrate a minimum in the energy (middle surface in Fig. [T]) and this makes it the 
yrast state, for this A and L. The optimal values of the sigmas are a = 0.72, <j\ = 0.70, the 
minimum energy per particle for these values of sigmas is eo = 1.37 and the s-depletion is 
d\ = 0.33. However the energy barrier, that prevents the system from collapse, is extremelly 
low, h ~ 10 -3 , making the state only marginally metastable. On the other hand, the I^/^q) 
excited state of the system exhibits a clear minimum (energy barrier height h = 0.23), with 
energy per particle eo = 1.60 and s-depletion di = 0.48 at the optimal values of the sigmas 
(T = 0.82, <Ti = 0.81 (upper surface in Fig. [Q). 

For L = a metastable fragmented state can decay by two channels. The first, as 
mentioned above by tunneling through the barrier. The second, by coupling to lower surfaces 
with the same L = which do not have a minimum. Since all these surfaces do not have a 
minimum are energetically far below, the coupling between them is not expected to induce a 
quick collapse. Consequently, metastable excited states with L = 0, with parameter values 
tuned at the collapse region of a GP state, exist at higher energies. 

B. Ground states for various angular momenta L 

Next, we perform the same analysis as in section UlI AI for the system of iV = 120 bosons, 
this time over states of significantly higher angular momentum. Precisely we choose states 
with L = 52, Ml = 0. We recall that the maximum allowed quantum number for the total 
angular momentum, within the present analysis, is L max = N = 120. We want to compare 
the stability and the properties of the two systems, namely that of L = to that of L = 52. 
The energy surface e(cr , <7i) as a function of the scaling parameters cr , &i is plotted in Fig. 
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[2j A clear minimum can be seen, at o"o = 0.82, o\ = 0.81 and eo = 1.59 manifesting a 
metastable ground state with L = 52, for the same system whose L = ground state is 
found to be collapsed. We should stress here that the state we examine is the lowest in 
energy state of this L and so this makes it the ground (yrast) state of the problem. 

In Fig. [3] we plot the energy surface of the same state, examined above, for different 
values of the interaction strengh, Ao = (—0.010, —0.056, —0.100). For small values of Ao = 
(—0.010, —0.056) the energy surface exhibits a clear minimum, with its energy barrier being 
higher than in the case of Ao = —0.0842. In the third picture, the energy surfaces shows 
no minimum, meaning that this state is collapsed, though the critical value Aq^, 52 is much 
higher than the corresponding Aq" of the L = state. 

Following Fig. [U a plot of energy surfaces of ground (yrast) states, i — 1, with different 
angular momentum and hence different stability behaviour, would be intuitive. If one would 
plot the energies of the group of ground states l^^rf 8 ) and l^^rf 8 ) on the (ao,^) 

plane, they would see that the resulting graph would look very much like that of Fig. [TJ 
This means that the energy surfaces of the pairs of states l^/^ 38 ) and l^/fsfo) as weu as 
I^^Ti 58 ) and |\l/^g) are almost the same, for all cxo,cxi. This coincidence is not an accident. 
Indeed, as we shall show later, one can find states that are very close — almost degenerate 
— in energy but have different angular momentum quantum number L (see discussion at 
the end of Sec IIIIC1D . 

As a direct generalization of the above, we can say that if, for some Ao the ground state 
with L = of the A-boson system is collapsed, then there will be a ground state with angular 
momentum L > 0, large enough so as to survive the collapse. Further, if the interaction 
strength is increased, past some new critical value, this state will also collapse. 

C. Analysis and structure of the energy surfaces 

To thoroughly analyze the properties of the MB states we examine the findings of the 
previous sections under the light of the natural orbital analysis and the use of RDMs. For 
given ground and excited metastable states we want to answer on: (i) what the natural 
occupations are, (ii) how much fragmented the states are and (Hi) how much they deviate 
from MF states, in a range of the parameters o~o,o~i as well as Xq,L,Mi. The systems 
examined in this section consist of N = 12 and N = 60 bosons and the interaction strength 
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is set to Ao = —0.842 and to Ao = —0.1684, respectively. 
1. Fragmentation 

As mentioned, due to the symmetry of the Hamiltonian, the natural orbitals of Eq. ( 122]) 
coincide with those defined in Eq. fl9]), for all Ao,o"o,o"i. It is interesting to see how the 
occupations pi, defined in Eq. ( 12~4"|) . of the ground and excited metastable states of definite 
L, vary in the (cr ,o"i) plane or change with A . Unlike p 2 , p 3 or p 4 , the quantity pi (or di) 
is invariant — for given L — under changes of the quantum number of the operator L z . 
Furthermore, as long as solely ground states are considered, i.e., i — 1, p\ determines the 
total angular momentum L. These properties make p\ a quite informative and representative 
quantity of the state 

For a system of N = 60 bosons in the ground metastable state with L = 26, = 2, we 
calculate the depletion of the 0i-orbital (s-depletion di), see Eq. ( )25l) . as a function of the 
parameters o"o,o"i. In Fig. [5] we plot the contour lines p\ = const., versus the parameters 
o"o,o"i. The energy landscape of this particular state, for this choice of parameters would 
look very much like the one of Fig. [2j To allow a monitoring of the energy surface, we also 
plot in Fig. [4] the contours (light grey) of constant energy e. The dashed line is the highest- 
in-energy contour that corresponds to a metastable state. It splits the graph into four parts; 
in the upper right one the 'trajectories' are bounded, while they are not in the other parts of 
the space (hyperbolic trajectories). Thus it resembles a separatrix of a phase space, whose 
trajectories meet asympotically only in a saddle point. The energy per particle has a local 
minimum €q = 1.55 at ao = 0.81, a± = 0.79 and at this point the s-depletion is found to be 
0.45, i.e., 55% of the particles of the system are excited to the orbitals <f)j, j = 2, 3, 4. 

Of special interest is also the change of the s-depletion as the system moves towards the 
collapse. To make this evident we have plotted on Fig. [4] an arrow marking the 'collapse 
path', i.e., the line that connects the minimum (green dot) with the saddle point (green 
square) of the energy surface, i.e., the maximum of the energy barrier. Along this path the 
system moves over the energy barrier towards collapse and it crosses contours of different p\, 
as collapse takes place the s-depletion of the state increases. We note that for large values 
of the scaling parameters, i.e., Uq,iji ^> 1 the s-depletion remains practically unchanged. 

Every state |\I/ L ) with definite angular momentum L is (2L + l)-fold energetically degen- 
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erate, due to the quantum number Ml- This means that the energy landscape of a state 
would not feel any change in Ml- Recall from Eq. ( 117)) that the eigenvalue of L z of a 
permanent |$) is simply Ml = n 4 — n 2 . Similarly it can be shown that, for a general state 
|\&), M L = Pi — P2 holds. Non surprisingly, this suggests that the occupations Pi,p3, i.e., 
the occupations of the two mi = orbitals, do not contribute to the z-projection of the 
total angular momentum L. However, as the quantum number Ml of a state with a given 
L varies, only the occupation p\ remains unchanged, while ps varies accordingly to keep the 
total number of particles fixed, i.e., p% = N — p 2 — Pi — Pi = N — Ml — 2p 4 — const.. This 
behaviour is depicted in Fig. [5] for a system of iV = 12 bosons in its ground state, for L = 6 
on the first and L = 11 on the second panel. On panel [S^a), at point Ml = 1, the first 
from above line (blue) corresponds to the occupation p 1; the second (yellow) to p 3 , the third 
(purple) to pi and the fourth one (green) to p 2 . On panel [5^6) the sequence is p3, pi, p 2 , pi, 
with the same coloring. The occupation numbers presented here are calculated at the op- 
timal values of 00,01, that minimize the total energy of the system. Both the energy and 
the optimal Cj, i — 0, 1, are invariant under changes of Ml- By comparing the two panels 
we see that the same pattern on the changes of the occupations is repeated, with pi fixed 
at different values; at L = 6, pi ~ 6 while at L = 11, p\ ~ 1. We infer that the behaviour 
of the occupations against Ml is a general feature, independent of L, N. 

In Fig. [61 for a system of N = 12, we show how the depletion di varies with increasing 
absolute value of interaction strength. In the left panel the dotted lines correspond to the 
six excited states, i = 2, . . . 7, of L = 0. The solid line marks the depletion of the lowest- 
in-energy metastable state (ground state), at each value of Ao, at the optimal ao.ai. The 
successive 'jumps' of this line take place at the critical values X l where the state collapses. 
Thus the plane of the figure is divided into the right 'collapsed half-plane' and the left 
'metastable half-plane'. Similar tendencies persist for states of different angular momenta. 
This is shown in the right panel, where we plot three curves that correspond to MB states 
of different angular momenta; the lowest one with L = 0, the middle one L = 2, and the 
upper one with L = 8, all with Ml = 0. Each curve is the value of the s-depletion of the 
lowest-in-energy metastable state with specific L against A = |Ao|(iV— 1). For low interaction 
strength (A < 7) the ground state is the state with L = and (almost) zero fragmentation. 
For larger values of the interaction strength the condensed state cannot support a metastable 
state anymore. Though, the first excited (and fragmented) state \^f= 2 ) is found to be non- 
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collapsed. 

An examination of the s-depletions of the different ground states of the right panel of 
Fig. fallows one a comparison of the respective energies; indeed, two states with the 
same s-depletion are expected to have the same energy. For example, the s-depletions of 
the states |^/^°) and l^f^f) (first and second from below lines, respectively) are very close 
to each other for the whole range of Ao that they exist and their energies -E^l 1 !/^ )] and 
^[l^^] 2 }] are found to behave accordingly. In fact, those two states belong to a family of 
states {\^i^ k )}k, whose members, defined by: 

i k + ^ = q, q eN*, (29) 

have, for A = 0, the same energy, i.e., 

E[\*t=£ k )] Xo =° const. (30) 

for all possible That is, all the states with L = L^,i = q — L/./2, for some positive 

q G N*, are degenerate in the absence of interaction. The degeneracy of such a group of states 



has been already noted in Ref. 48] and subsequent works. However the states considered 
there are those of Ml = L and hence the description becomes essentially two dimensional. 
In the case of A < and L\ < L 2 Eq. f )30|) transforms to: 



m^)] A <° (31) 

Namely, the decrease in the energy is larger in the state with the lowest angular momentum, 
when the attraction is switched on. This behaviour can be seen in the comparison of the 
states of different angular momentum L, on the right panel of Fig. |6j 

Summarizing, we see that the s-depletion is an informative quantity of the state, as it 
reveals information about the energy and the angular momentum, that carries. The 
s-depletion of a particular metastable state remains almost fixed for a\,ao ^> 1, while it 
changes rapidly as the system is driven to the collapse region of the surface. The s-depletion 
does not depend on the angular momentum Ml- Among states with different symmetries 
(quantum numbers) that are energetically degenerate at A = 0, the attractive interaction 
favours energetically the one with lower L, hence smaller M^-degeneracy. 
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2. Variance 



Besides the s-depletion of the condensate, the variances r« and |r|, defined in Eqs. (1271) and 
give information about both the structure of the stationary states and the dynamical 
behaviour of them. Although the calculation of time-dependent states are beyond the scope 
of this work, one can, based on the present results, comment on the expected dynamical 
stability of the states. In a fully variational time-dependent multi-configurational approach 
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52j both the permanents and the expansion coefficients are time-dependent, i.e., 



Y^u Ck(t)\&k(t)}- As shown in Ref. [39] the expansion coefficients Ck in |\&) = X/i^l^i) 
comprise a Gaussian distribution on their own, of width characterized by variance \t\. So, a 
state with a large value of |r| will include a large number of coefficients Ck in its expansion. 
For this reason it is expected to be dynamically more unstable than a state with small |t|. 

To study the variance of the states \^ L ) we plot in Fig. [7] the contours of fixed variance 
\t\ = const, on the (cro,o"i) plane for a system of N = 60 bosons in the ground state with 
L = 26, Ml = 2. We also draw the 'collapse path' (arrow) as defined before, the minimum 
(dot) and the saddle point (square) of the energy surface as well as the contours of constant 
energy e (grey lines). At the minimum of energy, at point e = 1.55, <t = 0.81, ui = 0.79, 
the variance of the system is \t\ = 4.79. As the systems moves along the 'collapse path' on 
the energy surface it crosses contours of different variance |r| towards larger values. Since a 
zero (or almost zero) value of |r| is indicative of a MF state, we see that the system moves, 
in this way, towards less and less MF states. On the other hand, for large values of the 
scaling parameters, i.e., o"i,o"o ^> 1, the variance \t\ remains practically unchanged. 

We have also examined the case of a non collapsed GP ground state of zero angular 
momentum. For values of parameters N = 12 and Ao = —0.5052 the ground state of the 
system is the condensed state with L = and the variance |r|, as well as the s-depletion 
di, at the optimal o"o,o"i is almost zero. The same as before scenario is found to hold; in 
a neighbourhood of the minimum of the energy, in the (o"o, Ci) plane, the variance remains 
very close to zero but as the system moves over the energy barrier the variance grows larger, 
i.e., the system moves towards non MF states. The same happens to the s-depletion d\. 
Note that, in all cases, the minimum value of \t\ and the optimum one (i.e., the value of \r\ 
at the minimum of energy) do not coincide. 

In Fig. [8] we plot the change in variance |r| of ground states \^ L ), against the quantum 
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number L, for six different values of the interaction strength Ao- The number of particles 
is iV = 60 and L varies from L = to L = 58 or < L/N < 0.97. As we increase 
the value of |Ao| the L-states, starting from L = upwards, collapse and hence cease to 
exist. We denote with L min the minimum value of L with which, at a given value of Ao, 
a metastable ground state of angular momentum L min can exist. For small values of |Ao|, 
where L TO j„ — 0, the variation of the states increases monotonously with L. For larger values 
of | Ao | (Ao < —0.15) a minimum in the curve r(L) appears, at a point L > L m i n > 0. The 
variances for all different values of the interaction strength meet at one point, as L — > N. 
Generally we detect two competing tendencies on \t\ as L increases; first, since the size of 
the configuration space N p drops linearly with L (N p = 1 when L = N) the number of 
coefficients in the expansion of Eq. ([5]) decreases with L and so will \t\. On the other hand, 
as L grows larger, the configurations $ include more basis-functions $ in their expansion 
and hence their variance |r|^ increases. The 'dominance' of the one or the other tendency 
seems to be conditioned by the value of the interaction strength Ao- However, for large 
values of L, the dependence of |r| on Ao is not significant. 

We, next, study the dependence of the variance |r| of the states |\I/ L ) on the quantum 
number Ml- We recall that the maximum angular momentum L max that a MB state can 
possess is, due to the orbital subspace used here, always equal to the total number of 
particles N. The (2L + 1) M^-states, of different z-projection of L, make every L-eigenstate 
(2L + l)-fold degenerate. 

In Fig. |9] we plot the variance |r| as a function of Ml for various states. For systems 
of (a) N = 12 and (6) = 60 bosons we choose three different ground states l^^) with 
L = 5,8,11 and L = 26,40,58 (first and second panels, respectively). In the figure, at 
Ml = 4 for the left and Ml = 20 for the right panel, the lowest, middle and upper curves 
correspond to the lowest, middle and upper values of L, respectively (blue, purple and 
yellow colors). As the quantum number Ml increases the variance \t\ drops, contrary to the 
fact that the size of the configuration space N p does not depend on Ml (see Appendix IA1 . 
However, the size of the expansion of the basis functions $ scales like (N — \Ml\) 2 and this 
results in the decrease of the variance |r |^ of each of the functions as Ml increases. In the 
'edge' of each L-block, where Ml = ±L, ±(L — 1), the variance takes always its minimum 
value (see also Appendix IB 2j) . If, further, L = N and Ml = ±L = ±A^ the variance |r| 
is zero, since there is only the permanent |0, N, 0, 0) (or |0, 0, 0, N)) that contributes to the 
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state |*). 

Note that the shown dependence of the variance \t\ on Ml is connected to the size of the 
(truncated) space of one-particle basis functions that we use. In similar calculations over 
an extended (i.e., less truncated) </>-space, there would be more terms in the expansions of 
$ and the variances shifted to higher values. However the general tendencies, as shown in 
Figs. |8] and [9] are not expected to change. 

In this section we have studied the dependence of the variance \t\ of a state on the 
parameters o"o, <J\ and the quantum numbers L and Ml- Generally, as the system moves 
towards collapse (i.e., <Ji,a — > 0) the variance |r| increases. Moreover, the variance as a 
function of L can increase monotonously or exhibit a minimum, depending on the value of 
Ao- The variance |r| decreases with increasing Ml- 

IV. ANGULAR MOMENTUM AND COLLAPSE: MANY-BODY VS. MEAN- 
FIELD 

As already discussed, any three dimensional attractive condensate is expected to col- 
lapse when the product A = |A |(iV — 1) exceeds a critical value A cr . However, fragmented 
metastable states can survive the collapse for a much greater value A > A cr . In this section 
we examine the behaviour of MB states \^ L ), as well as these of the MF states |$) of various 
angular momenta — exact or expectation values — in the onset of collapse. Combining the 
findings of the previous discussion we show how the angular momentum can stabilize an 
overcritical condensate. We first discuss the impact of angular momentum on the stability 
of MB states. We then give an account of the estimated angular momentum within the MF 
approximation by deriving relevant quantities (expectation value of the angular momentum 
operator) that will allow us comparisons with the MB results. 

A. Many-Body predictions 

In the previous section, we described the structure of MB states that have a definite 
angular momentum ^ L ^ N. We showed that, generally, these states are fragmented 
and, moreover, are non MF states. This suggests that a MB state |\1/ L ) with definite L can, 
depending on its s-depletion and the value of |Aq|, survive the collapse. Additionally, the 
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condition [H, L] — necessitates the conservation of the total angular momentum and thus 
the stability of the state \^ L ). 

Figure [10] summarizes and aggregates the main results of this work. We first focus on the 
upper connected dotted lines, which are the results for the MB states. For systems of different 
particle numbers N = 12, 20, 60 and 120 (see the legend of the figure for the correspondence 
to the different colors) we plot the s-depletion d\ versus the quantity A = |Ao|(iV — 1). Each 
plotted point, at each value of A, is the depletion d x of the ground (yrast) state 1^^) of some 
angular momentum L which is still non-collapsed. As the absolute value of the interaction 
strength increases, the lowest-in-energy states \^ L ) start to collapse. The energies and 
occupations (depletions) are calculated at the optimal values of the parameters 00,01. As 
we have already seen in Sec. IIII C 1[ at a given A , the s-depletion of a MB ground state 
gives also the angular momentum L/N of this state. Qualitatively, for the ground state of 
each L-block, one can write 

^ = 1 - % + O (r(A )) = d x + O (r(A )) , (32) 

i.e., the angular momentum of a ground state and the depletion of it differ only to some 
term 0(t), that depends on the fluctuation (variance) of that state, which in turn depends 
on the strength of the interaction. In a non-interacting system the fluctuations are zero and 
d\ = 4 exactly. 

Interpreting the results of Fig. [10] we can say that for any value of the factor A there 
will be some L > such that the (ground) state l^^i) is metastable. The critical angular 
momentum L increases monotonously with A. The stability behaviour seems not to depend 
significantly on the number of bosons in the following sense: for small particle numbers the 
curves of Fig. [10] are slightly different, while for A^ 3> 1 all curves converge, rendering in 
such a way the obtained results universal and independent of a particular choice of Ao or N. 

B. Mean-Field predictions 

Any MB state as we saw, is an eigenfunction of the operator L 2 . At the MF level 

however every state |$) of the system is represented by only one permanent, Eq. (J2]). Hence, 
with the exceptions of states with Ml = ±L, ±(L — 1), a MF state |<&) is by construction 
incapable of describing eigenstates of L 2 (see Appendix IB 21 for the possible MF states that 
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are eigenstates of the total angular momentum operator). This incapability comprises a 



major difference between the two descriptions. Within the multi-orbital BMF 



2l| theory 



the occupations n, of each orbital of the ground state are varied to extremize the energy 



functional of this state. However, in the description of excited states [24( they serve as 
parameters that are externally determined. In such a way one is free to choose the values for 
the set of the occupations {ni, n^, n 3 , n 4 } or {ri2, n 3 } for given depletion dn and total particle 
number N. So, for example, the choice ri2 = ns = ^ 0, made in Ref. 24j, guarantees the 
sphericallity of the one-particle density [i.e., p(r) = p{r)\, but breaks the L-symmetry of the 
state. We recall that in the present MB approach the natural occupation numbers, for all 
the ground and excited states, are determined variationally from the eigenvectors C of the 
optimized Hamiltonian matrix H, see Eq. f[2~41) . As a result, the rotational symmetries of 
the system are restored. 

So, what is the angular momentum that MF states have? It is a matter of fact that 
at a MF level one can only speak of expectation values and not exact values/quantum 
numbers of L. It can be shown (Appendix IB 21) that the expectation value (L) of the 
angular momentum of a MF state with equally distributed excited bosons ri2 = n 3 = is 
the statistical average (mean) of the exact total angular momentum of the MB states with 
the same value of depletion d±: 

Lmf = (Lmb)^, (33) 
where L MF (L MF + 1) = (L 2 ). So, in accordance to its name, the mean-field state can 
provide only the mean angular momentum of the corresponding (i.e., same d±) MB states. 
Furthermore, one can calculate (Appendix IB 21) the average momentum over, first, all the 
MF states |<&) and, second, over all MB states. It turns out that they are connected through: 

N 3 

(Lmf) all states = = ^J^(Lmb) all states- (34) 

So, the average angular momentum over MF states is found to be ^ ~ 1.22 times higher 
than the average one over MB states. 

Since, within the MF theory, the states |<3>) do not possess a definite quantum number 
L we cannot write any exact correspondence between the depletion di of the state and the 
angular momentum L, as we did in the case of MB ground states. Instead, we can use Lmf 
and relate it to d\ through: 

X ~rVr =o^l- (35) 
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This result is taken in the limit N ^> 1 (see Appendix IB 2j) . Note that it does not depend 



on the value of Ao- This reflects the absence of fluctuations on a MF state, which do depend 
on the interaction strength Ao- 

How is Lmf related to the stability of the condensate? Recall first, that a system will 
survive the collapse if, for a given Aq, the number of particles that occupy the s-orbital stays 



below a critical number N cr 



60j . This however does not forbid the total number iV of bosons 
of the system to be larger than this critical number. Indeed a system can exist in a state 
with ri\ < bosons occupying the s-orbital and N — n\ occupying higher-in-energy or- 
bitals. More precisely, any excitations of bosons to p-orbitals may increase the total energy 
of the system but will contribute to the total stability of it, since the excited p-bosons 'feel' 
less the interaction energy than the s-bosons. This is the reasoning behind the metasta- 
bility of fragmented states with an overcritical number of bosons, already demonstrated in 
Ref. 24|. Here we further show that a MF state |$) with non-zero expectation value of 
angular momentum Lmf > exhibits fragmentation, Eq. ( 1351) . which increases the overall 
stabilility of the system. However, the impact of the angular momentum in the stability 
of the condensate is overestimated at the MF level. A comparison of Eq. ( 135]) with the 
corresponding MB one, Eq. ( 132]) . convinces us of this claim. 

To allow a better comparison to the MB results of Fig. [10] we plot on the same graph 
the data obtained from the MF states (second group of dotted unconnected lines in Fig. 
[LO]) . More precisely, for systems of N = 12,20 and A^ = 60 bosons, we plot at each value 
of |Ao|(A/" — 1) the s-depletion d± = 1 — pi/N, with p% now given by the critical number 
of particles N cr (i.e., maximum number of particles so that the state |<E>) is not collapsed) 
calculated from the relation: 

(-64 ■ 5 3 / 4 - 128 ■ 5 3 / 4 f + 300Aq - 375A f ) 
cr cr [_ 4 + 13(f) 2 ] (16 -5 3 / 4 - 75A ) ' 1 J 

at a = a Xl with A = ^ (f ) 1/2 f , Ng p = 1 - (±) 1/4 and n = n 2 = n 3 = n A the 
occupations of the p-orbitals. Equation f ]36]) in the limit A^ ^> 1 gives back Eq. (7) of Ref. 



24| . Also, the numerical MF calculations for the critical numbers of a A^ = 120 bosons 
system, without using the assumption cr = &i, are presented in Fig. CG)] with the 'boxed' 
line (blue). The second from below continuous line (dark yellow) determines the angular 
momentum expectation value L MF /N [Eqs. (135]) and (1B29[) ] . over MF states. The lowest 
continuous line (red) of Fig. [10] is the calculations from the Gross-Pitaevskii theory. In this 
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case one has to identify 1 — with 1 — where Nf r p is the maximum number of bosons 
that, for a given Ao, can be loaded in a GP state without collapse. Here we use N = 60 
bosons. The total particle number N^ r p is, of course, the number of s-bosons of the system. 
Obviously this critical number is decreased, as we move to the right of the x-axis of the 
diagram and hence we call this curve the 'critical GP'. 

The 'bands' of MF and MB states depicted in Fig. [10] substantially deviate one from 
each other at small and moderately larger values of A. This is nicely manifested in the 
difference between the MF and MB predictions of the collapse of the L = ground state. 
The collapse of the MB state appears to happen at a smaller value of the product A than 
the one that the MF theory estimates. This reflects the overestimation of the impact of the 
angular momentum within the MF and puts the MB prediction closer to the experimentally 
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measured values of A (see Ref. js| and also the discussion in Refs. 
discrepancies between MF predictions and experimental values of the critical numbers and 
the collapse times). 

We see that the form of the curves for the s-depletion of the MF states seems not to 
be affected from the number N of the particles of the system. The various plotted MF 
curves for different N, like the MB ones, tend to converge for N ^> 1, making thus the 
described stability behaviour a universal and independent of iV phenomenon. For the MF 
case convergence has been noticed already for N ~ 10 2 bosons. Note, though that, unlike the 
MB states, the MF ones with L = collapse all at the same critical value A cr = |Ao, C r|(^ — 1), 
regardless of the total number N of bosons. We also see in Fig. UHla divergence of the angular 
momentum Lmf (dark yellow line) from the s-depletion of the MF states (dotted lines); this 
is exactly the relation of the two quantities, that Eq. f )35|) provides. The 'critical GP' curve 
significantly diverges from both the multi-orbital MF and the MB predictions. 

Conclusively, we presented a way, Eq. ( 133]) . to connect the angular momenta of a MF 
state of the form |m, n, n, n) to that of the MB states with the same depletion d\. A non- 
zero angular momentum will result in a fragmented condensate [Eq. ( |35]) ] which in turn 
will render the system more stable, with respect to the parameter A. Those results are in 
agreement with the MB ones of the previous section. 
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V. SUMMARY AND CONCLUSIONS 



In this work we constructed many-body states with definite angular momentum quan- 
tum numbers L and M L , for systems of N isotropically trapped bosons in three dimensions, 
interacting via an attractive two-body potential. These many-body states are written as an 
expansion (configuration interaction expansion) over orthogonal many-body basis functions 
(permanents). We represented the Hamiltonian and angular momentum operators as ma- 
trices on this basis and we looked for the states that simultaneously diagonalize them. In 
this representation the Hamiltonian has a block-diagonal form, with each block consisting 
of many-body states, with the same eigenenvalue of angular momentum. The one-body 
basis functions that we used are the wave functions (s- and p-orbitals) that solve exactly the 
linear (non-interacting) problem, each scaled under a parameter <jj, which we determined 
variationally. The rotational symmetries as well as symmetries under spatial inversion that 
the one-body basis functions possess are also present in the many-body states and reduce 
significantly the size of the configuration space. Due to the truncated one-particle basis set, 
the total angular momentum is restricted to < L ^ N. To our knowledge this is the first 
time that a fully three-dimensional Bose gas in isotropic trapping potential is studied, with 
the many-body wave function of the system expressly written as an eigenfunction of both 
total angular momentum operators L 2 and L z , for L > 0. 

For a value of the parameter A = |A |(iV — 1) such that the L = ground state of 
the system is collapsed, we have plotted the energy per particle e(<7o,<7i) of the ground 
and the excited many-body states, as a function of the parameters a ,ai. We have shown 
that metastable excited states of the same angular momentum can exist. Furthermore, for 
the same system, we demostrated the existence of metastable ground states with angular 
momentum L > that can survive the collapse. These states would also collapse, if the 
(absolute value of the) interaction strength is further increased. The examination of the 
above states, in terms of the natural orbital analysis, revealed that the states are fragmented, 
with a substantial number of particles being excited to the p-orbitals. 

We discussed why the s-depletion of a many-body state bears information about the 
energy and the angular momentum of |\&). We found that the s-depletion of a metastable 
state remains practically fixed for <7o, o\ ^> 1, while it changes rapidly as the system is driven 
to collapse. We have shown also that the z-projection of the angular momentum M L does 
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not affect the occupation of the first natural orbital. 

We have studied the dependence of the variance |r| of a state on the parameters <xo,o"i 
and the quantum numbers L and M L . We saw that along the 'collapse path' the variance 
increases. The variance as a function of L — depending on the value of Ao — can increase 
monotonously or exhibit a minimum. We also found that as the quantum number 
increases the variance |r| decreases. 

To further investigate the impact of the angular momentum on the stability of the system, 
we plotted the critical s-depletion of the metastable ground states \^f =1 ) (yrast states) as 
a function of the quantity A = |A |(iV — 1). We showed the connection of the s-depletion 
to the critical angular momentum L, in both the mean-field and the many-body cases. We 
demonstrated that for any value of the factor A there is some angular momentum L > such 
that the (ground) state \^i =1 ) is metastable. The critical angular momentum L increases 
monotonously with A and this behaviour is found to be independent of the particle number 
N, as long as JV 3> 1. We derived analytical relations for the expectation value of the 
angular momentum of a mean-field state, with equally distributed excited bosons, which 
allowed us to compare it with the corresponding results from the many-body approach. We 
have further demonstrated that the angular momentum of this mean-field state equals the 
average angular momentum of many-body states, with the same s-depletion. 

Conclusively, we can say that for any particle number N and interaction strength A of 
an attractive condensate, there is some well defined quantum number L of the many-body 
angular momentum operator L 2 such that the ground state of this system is metastable, i.e., 
exhibits a clear minimum in energy as a function of the shapes of the orbitals. Moreover, 
since the total angular momentum of the system is conserved, once the system is prepared 
in a ground state with L > it can survive the collapse, and that for a particle num- 
ber/interaction strength much beyond the corresponding ones of the L = ground state. 
We hope that our results will stimulate experimental research. 
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Appendix A: Size of Fock Space 



The total number of the iV-body basis functions (permanents) that can be constructed 
over a basis of M one-particle wave functions of Eq. is 




N - \ - I - ( M + N -^ (A1) 

I - I " N\(M-1)\ ' ( } 

which for M = 4 becomes 

1 N 3 
N P = -(N+1)(N + 2)(N + 3) ~ — . (A2) 
o o 

Using the symmetries of the system we can reduce significantly this number and hence the 
complexity of the problem. Without loss of generality we assume that the particle number 
N and the quantum number Ml are even integers. 

a. Total angular momentum L z : Since [L z , H] — the state |tf ) = J2i dl^i) 

can be 

chosen to be a common eigenfunction of the two operators. This transforms H to a block 
diagonal form, with every block consisting of states of distinct Ml- The number of states 
in a block with some Ml is 

N P = (?L±id^ < ^. (A3) 

b. Parity fl^r) = \&(— r): Similarly, [n, H] = and T-L block diagonalizes into two 
blocks, each with distinct parity II = +1 or II = — 1. The number of states |^) in the block 
with II = 1 is 

[N + 4 + 2M L - 3M L H(M L )} [N + 2- M l H(M l )} < 

P g ~ g ' v / 

where H{x) is the unit-step function. 

c. Total angular momentum L 2 : Last, the commutator [L 2 ,H] = 0, diagonalizes the 
matrix % into blocks of states that have definite angular momentum quantum number L. 
The number of states |\I/ L ) in the block with some L is 

N-L+2 N 
~2 ~ ~2 

where L is the quantum number of L 2 and here it is assumed to be an even number. In case 
L is odd Eq. f ]A5j) should read: N p — (N — L + l)/2. Note that these relations hold for any 
M L . 
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N P = ^— < -, (A5) 



Appendix B: Angular momentum in Many-Body and Mean-Field theories 



1. Many-body eigenstates of the total angular momentum operator I? 

We now return to the question of explicitly finding the eigenstates of the operator L 2 , as 
discussed in Sec. Ill Bl 

A general permanent |$) = |n), representing a system of a total number of bosons 
N = rii + n 2 + n 3 + n 4 and z-projection of the angular momentum Ml = n 2 — n 4 , takes on 
the form: 

\n) = \ni,n 2 ,n 3 ,n 4 } = \N - 2n 2 - n 3 + M L ,n 2 ,n 3 ,n 2 - M L ), (Bl) 

where n 2 ,n 3 are integers, such that M L H(M L ) ^ n 2 < (N+M L )/2, ^ n 3 < N-2n 2 + M L , 
where H(x) is the unit-step function. An expansion over these (orthogonal) permanents 
|$) is: 

l*>= J2 C ^ ,1*0 (B2) 

n 2 ,n 3 

where n 2 , n 3 run over all possible permanents of Eq. (IB II) . Acting operator Eq. ( fl3l) on Eq. 
( 1B2!) we get: 

L 2 |v[/> = L 2 Cn 2 , m \n) = AJ2 C n 2 ,n s \n), (B3) 

712,713 n2,ns 

or 

A J] C n2tn:i \n) = C n2>na (A{n 2 ,n 3 )\n) + B{n 2 ,n 3 )\n + 2) + T(n 2 ,n 3 )\n - 2)^, (B4) 

where A = L(L + 1) are the eigenvalues of L 2 , \n + 2) = |ni,n 2 — l,n 3 + 2,114 — 1) and 
|n — 2) = |ni, n 2 + 1, n 3 — 2, 714 + 1), i.e., they are the double 'excitations' of the permanent 
\n). The functions A,B,T are: 

A{n 2 , n 3 ) = n 2 (n 3 + 1) + n 3 (n A + 1) + n 3 (n 2 + 1) + n 4 (n 3 + 1) + (n 2 - n 4 ) 2 , 

5(n 2 , n 3 ) = 2 [n 2 n 4 (n 3 + l)(n 3 + 2)] 1/2 , (B5) 

r(n 2 , n 3 ) = 2 [n 3 (n 3 - l)(n 2 + l)(n 4 + 1)] 1/2 . 

The problem is focused in calculating the coefficients C mjTl3 such that Eq. f lB3j) is fulfilled. 
We will show how one can reduce this equation to a simpler form. By multiplying Eq. (IB4j) 
with (n| and using orthogonality of permanents and the relation 

T(n 2 - k, n 3 + 2Jfe) = 5(n 2 - Jfe + 1, n 3 + 2A; - 2), k G N, (B6) 
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we obtain: 

ACn 2 ,n 3 — A( n 2, n 3 )C n2tn3 + T(n 2 , n 3 )C n2+ i t n 3 ^ 2 + B(n 2 , n 3 )C m _i in3+2 . (B7) 

This is a homogeneous, second order recurrence (or difference) equation of the two indepen- 
dent variables 112,11,3, with known non-constant coefficients. 

In the above equations there are two free parameters n 2 , n 3 which are varied independently 
and these are also the independent variables of Eq. ( IB7j) . To reduce the dimensionality of 
the problem one can proceed by switching the representation of the permanents and their 
coefficients. Precisely, we can use a simpler representation for indexing the vectors |$) in 
the expansion of \^). Noticing that the action of the operator L 2 on a state of Eq. (1B2[) 
involves only permanents of the form 

\fk) = I Tii, a — i, (3 + 2i, a — i + Ml), (B8) 

where a, [3 G N and — (N + Ml)/2 ^ i ^ —MlH(Ml), we can work with permanents of 
the above type only, for fixed a, (3. In fact the action of L 2 partitions the configuration 
space into invariant subspaces, with permanents of the form of Eq. ( 1B8I) . Permanents with 
a 7^ a' or (3 ^ (3' will not contribute to the same eigenstate This allows us to move 
from the two-parametric {n 2 ,n 3 }~ to the one-parametric {i}— representation. We write 
now again Eqs. flB2l) - flB7l) in the new representation. 
A general state becomes: 

\*) = Y,Ci\fli), (B9) 



where i runs again over all permanents (IB8I) . Similarly, acting operator Eq. ffT3l) on Eq. 
(1BH we get: 

L 2 \^) = L 2 Y,Ci\n l ) = A^QItO (BIO) 

i i 

= ^2Ci(Ai\ni) + Bi\ni + 2) + - 2)), (Bll) 



with Aj = A(a — i, f3 + 2i), Bi = B(a — i, (3 + 2i) and Tj = T(a — 1,(3 + 2%). Equations ( IB6 
and ( 1B7I) become: 



and 



^i+2k — B i+2 (k-i), (B12) 

(A { - A)Ct + r.CU + B t C l+l = 0, (B13) 
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respectively. The above is a homogeneous second-order recurrence (difference) equation of 
one independent variable [cf. Eq. (IB 71) ]. 

For some choices of the parameters a, (3, A, M L Eq. (1B13P can be easily solved. In partic- 
ular, for a = 0, (3 = N and Ml = we obtain: 

I -i\2N+l~ir,N+2-2i-p( N , 9 ,^ r ; N+l \ 

A = (L = 0), C t = tll r{ J/ 2 J + ^^ Co, (B14) 

V ij, (6i-JV-6)r(i-i/2) °o, 

a _«(/-_ z\ r - (-i) 2JV+1 -'2^+ 5 - 2 '(jV-2)jvr(f+ 2 - t )r(iv±i) , , 

l\ u yu 6), (140i 2 -60(Ar+5)j+3Af(Af+18)+160)r(i-l/2) ° ' V^ 1 ^) 

where Cq is to be determined from the normalization condition ^\ |Cj| 2 = 1 and —N/2 ^ 
i ^ 0. Equations ( 1B14l) -( lBT6l) give three of the states $ L , that are eigenstates of L 2 and 
belong to the rotated basis {<J>f } of Sec. Ill Bt 



2. Mean-field and average many-body angular momentum 

We show here that the total angular momentum of a mean-field state, with equally 
distributed excited bosons n 2 = = n 4 = n is the statistical average of the exact total 
angular momentum of the many-body states with the same depletion d\ — 1 — i.e., 

Lmf = (L MB ) dl . (B17) 

Recall that ^ L ^ N. Then, for a total number of iV bosons there are N + 1 blocks 
(L-blocks) of the Hamiltonian matrix H, each with a distinct value of L. We want to 
calculate the average angular momentum (L), among states \^ L ) ni with a given natural 
occupation p\ = n\. We assume that in each L-block this occupation n\, as we move from 
the highest-excited state to the ground state, increases like 

nUk) = 2k, if L is even 

(B18) 

ni(k) = 2k + 1, if L is odd, 

where k e N indexes the state \^i =k ). The above relations hold exactly in the absence of 
interaction, i.e., Ao = 0, and in a satisfactory approximation when Ao ^ 0. Then, as we have 
numerically verified, each L-block with L < iV — n± contains exactly one state \^ L ) ni with 
the desired n\ (or very close to it). Recall that the size of an L-block drops linearly with 
L, as in Eq. (IA5j) . So there are N + 1 — rii L-blocks that contain one state \^ L ) ni . The 
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occupation n\, in the case of Ao = 0, is even in half of the blocks, odd in the other half ones. 
The total number of states with occupation n\ is: 

jv+i-m,2 

N(n x )= Yl + ( B19 ) 

i=0 

due to the L z degeneracy. To include only even (or approximately even) occupations n\ we 
sum on a step of two (the added term ',2' in the upper limit of the sum denotes that step). 
These states have total angular momentum: 

iV+l-711,2 

Ltotal= ( 2l + 1 ) L i- ( B2 °) 

i=0 

The quantum number Lj of each block simply increases like Lj = % and hence 

EM2i + l) N'(4N' + 7) 2 
= £(2i + 1) = W + 1) 3 (iV ^ (B21) 

where N' = N — n\. In the limit iV > 1 we get: 

(L MB ) dl ~ HjVdx. (B22) 
The average over all MB states, of all n\ simply gives: 

(L MB ) all states = J- (B23) 

On the other hand, a MF state, with equidistributed excited bosons: 

|$) = \ni,n,n,n), (B24) 

with ni + 3n = iV, has no well-defined angular momentum quantum number L (except from 
the single case \N, 0, 0, 0)). We can, though, calculate the expectation value on a state |<&) 
from Eq. (fl6l) as: 

(L 2 ) = (n|L 2 |ra) = ra 2 (ra 3 + 1) + n 4 (n 3 + 1) + n 3 (n 4 + 1) + n 3 (n 2 + 1) + (n 2 - n 4 ) 2 (B25) 
for a general permanent 

|$) = |7ii,n2,n 3 ,n4>, (B26) 

or 

(L 2 ) = An(n + 1) (B27) 
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for the permanent of Eq. (IB24I) . For comparison purposes, we define a pseudo-quantum 
number Lmf, such that 

L MF (L MF + 1) = (L 2 ). (B28) 

Hence: 



Lmf = - (-1 + Vl6n 2 + 16n+ l) = ^ 
For N 1, we get: 



1 fi 1 



(B29) 



Lmf = ~{N- ni )~ (L MB ) ni . (B30) 

So, indeed the angular momentum of the MF state (1B24I) equals, under the assumption 
N ^> 1, the mean angular-momentum of the MB states with the same s-depletion. Equation 
( 1B30I) immediately gives back Eq. (|35l) : 



iMF ?fl-=n=?4. (B3D 



AT 3 V N J 3 
Now, the average angular momentum over the permanents of Eq. (1B26|) with the same 
rii is ; 

<W>„, = -\ + /^-^iiV, (B32) 
for iV ^> 1 and AT > 3n 1; whereas the average over all the permanents of Eq. (IB26I) reads: 

N 3 

{Lmf) all states = = \Lmb) all states-, (B33) 

also at N > 1. 

Last, we prove the condition for a MF state of Eq. (IB26[) to be eigenstate of the an- 
gular momentum operator L , already given in Sec. IIVBI Let |$ L ) be a single-permanent 
eigenstate of L 2 of Eq. (TT3"j) . Then it must 

L 2 |$ L ) = L(L + 1)|$ L ), (B34) 

where L{L + 1) is the eigenvalue of L 2 for this permanent. Then from Eq. fll6p we get that 
the conditions: 

«3 = or n% = 1 and 

(B35) 

n 2 = or n 4 = 
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must hold simultaneously. From here it turns out that the permanents that can satisfy Eq. 
flB34j) are: 



|$ L ) = \N + M L , 0, 0, -M L ), with M L = -L, (B36) 

|$ L ) = \N- M L , M l , 0, 0), with M L = L, (B37) 

= \N + M L -l, 0, 1, -M L ), with M L = -L + 1, L > 1, (B38) 

= |jv - M L - 1, M L , 1, 0), with M L = L - 1, L > 1, (B39) 

where iV the total number of particles and Mx = n2 — the quantum number of L z , as 
usual. Thus we see that the only permanents that can be eigenfunctions of the operator L 2 
are the permanents with quantum numbers restricted to: 

M L = ±L,±(L-1). (B40) 

Unless A = 0, Eq. f]B40|) serves as a necessary but not sufficient condition, for a MF state 
to be eigenstate of both the angular momentum operators L 2 , L z and also the Hamiltonian 
H. In the case of Ao = there are no couplings among states with the same L and Ml and 
condition ( 1B40j) . hence, suffices to determine a MF eigenstate of L 2 , L z and H. The same is 
expected to happen for small values of Aq. 
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FIG. 1: (Color online) Energy landscape e(o"o>°"i) f° r a system of N = 120 bosons, Ao = —0.0842. 
Shown are energy surfaces of the ground (i = 1), the i = 20 and the i = 30 excited MB states, all 
with L = Ml = 0. The lowest surface corresponds to the ground state of the coherent system, where 
almost all 120 bosons reside in the s-orbital. It exhibits no minimum and hence the system collapses. 
The middle surface barely exhibits a minimum, at gq = 0.72, o\ = 0.70, with barrier height 
h = 3.83 • 10~ 3 . The natural occupations are p\ = 80.82, pi = 13.06, = 13.06 and p^ = 13.06. In 
the third surface a clear minimum in the energy, eo = 1.60, is shown, at cjq = 0.82, o\ = 0.81, with 
barrier height h = 0.23. The occupation numbers, at this point, are p\ = 62.13, p2 = 19.29, ps = 
19.29, pi = 19.29. All quantities are dimensionless. 
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FIG. 2: (Color online) For the same system as in Fig. [TJ i.e., N = 120, Ao = —0.0842 we plot the 
energy surface of the ground state with angular momentum L = 52, Ml = 0. Clearly there is a 
minimum in the surface, which manifests metastability of the system. Contrarily, when L = (Fig. 
[1]) the ground state is found to be collapsed. The minimum energy per particle is eo = 1.59, at point 
ao = 0.82, cxi = 0.81 and the occupation numbers p\ = 62.30, p2 = 14.30, p3 = 29.10, p4 = 14.30. 
All quantities are dimensionless. 
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FIG. 3: (Color online) Energy surfaces for the state of Fig. [2]for values of interaction strength Ao = 
—0.010, Ao = —0.056 and Ao = —0.100. The first two surfaces exhibit minima, i.e., metastability, 
while the third one does not and hence the system collapses. See text for more details. All 
quantities are dimensionless. 
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FIG. 4: (Color online) Change in the s-depletion of the condensate (d\ = 1 — pi/N) on the (<7o, 01) 
plane for a metastable ground state with L = 26, Ml = 2 of an N = 60, Ao = —0.1684 system. 
Plotted are contour lines of fixed p\ = const.. The minimum of energy (green dot in the plot) is 
eo = 1.55 at oq = 0.81, a\ = 0.79 and the s-depletion is d\ = 0.45. The saddle point (green square) 
on the energy surface gives the maximum energy that a metastable state can have. Along the 
'collapse path' (arrow) d± increases; by moving the system — over the energy barrier — it becomes 
more and more fragmented. All quantities are dimensionless. 
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FIG. 5: (Color online) Occupations with respect to Ml, for a system of N = 12 bosons, with 
interaction strength set to Ao = —0.842, in the ground state of angular momentum (a) L = 6 
and (6) L = 11. The occupations of the three excited natural orbitals differ significanlty when, 
inside the L = const, subspace, we increase the projection Ml of the angular momentum. The 
occupation numbers presented here are calculated at the optimal values of <jq, o\. Both the energy 
and the optimal dj, i = 0,1, are invariant under changes of Ml- Note that, for different values of 
L, the same pattern on the occupations pi persists, though p\ is fixed at different values, according 
to pi ~ N — L [see Eq. ()32jl ] . All quantities are dimensionless. 
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FIG. 6: (Color online) s-depletion, d\ = 1 — varying with the (absolute value of the) interaction 
strength for a system of N = 12 bosons. In panel (a) the dotted lines correspond to the ground 
and the six excited states of the L = block of the Hamiltonian. The solid line marks the lowest- 
in-energy metastable (yrast) state that was found at each Ao point. In panel (6), shown are three 
curves corresponding to the MB (yrast) states of different angular momentum L; the lowest one 
with L = 0,Ml = (blue), the middle one with L = 2, Ml = (magenta), and the upper one 
with L = 8, Ml = (yellow line). For weak interaction strength the ground state is the |^^°) 
state with (almost) zero fragmentation. For A = |Ao|(iV — 1) ~ 8 the lowest-in-energy state that 
survives the collapse is the fragmented state l^^ ) with d\ ~ 0.2. Its energy is very close to that 
of the ground state l^^ 2 ) of the L = 2 block. Compare also the three states l^^ 8 ), and 
at point A ~ 15 (see text for more details). All quantities are dimensionless. 
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FIG. 7: (Color online) Change in variance \t\ on the (o~o,o~i) plane for a N = 60 bosons system, 
for the same state as in Fig. SJ Plotted are the contours of |r| = const, as well as the contours of 
constant energy (grey curves). At the minimum of energy (green dot) eo = 1.55, <jq = 0.81, o\ = 
0.79, the variance of the system is |r| = 4.79. The arrow joins the minimum (dot) and the saddle 
point (square) of the energy surface, i.e., it depicts the 'collapse path'. As the system moves along 
this 'collapse path' the change of the variance grows moderately large (see text for more details). 
All quantities are dimensionless. 
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FIG. 8: (Color online) Change in the variance |r| of ground states l^^li) against the quantum 
number L, for six different values of the interaction strength Ao- The number of particles is N = 60 
and the maximum angular momentum is L max = N (in the diagrams up to L = 58). The curves 
shown are for metastable states which do exist. As the value of Ao grows larger the L-states, 
starting from L = upwards, collapse and hence cease to exist. L m i n is the minimum value of L 
that, at each value of |Ao| a metastable ground state of angular momentum L m i n exists. As an 
example L m i n is indicated by an arrow for the Ao = —0.16 curve. For small values of Ao, where 
L m in = 0, the variance of the states increases monotonously with L. For larger values of Ao a 
minimum of \t\ appears at some L > L m i n . All quantities are dimensionless. 



44 




FIG. 9: (Color online) Change in the variance |r| of MB states I^aLi ) against the quantum number 
Ml- Precisely, for a system (a) of N = 12 bosons and interaction strength Ao = —0.842 and (b) 
of ./V = 60 bosons and interaction strength Ao = —0.1684 we pick three ground states l^l^i) with 
different values of L. On the left panel, at Ml = 4, the bottom line (blue) corresponds to L = 5, 
the medium one (purple) to L = 8 and the upper (yellow) one to L = 11. On the right panel, 
at Ml = 20, the bottom line (blue) corresponds to L = 20, the medium one (purple) to L = 40 
and the upper (yellow) one to L = 58. In the 'edge' of each L-block of the Hamiltonian matrix H, 
where Ml = ±L,±(L — 1), the variance takes always its minimum value. For Ml = ±L = ±N 
the variance |r| is zero and the state is a MF state (see Appendix IB 2|h The shown decrease of 
the variance is attributed to the decreasing number of available permanents <3? that comprise the 
basis functions $ = as the number Ml increases (see text for more details). All quantities are 
dimensionless. 
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FIG. 10: (Color online) Stability plot for systems of different number of bosons N = 12, 20, 60, 120. 
As the interaction strength |Ao| increases the lowest-in-energy states start to collapse; here we plot, 
at each value of Ao, the depletion of the lowest-in-energy state, that is still non-collapsed. The 
connected dotted lines (upper 'band' of curves) are the many-body calculations. Every point of 
the MB plots corresponds to a ground (yrast) state l^^i) which, unlike the mean-field states, 
have a definite value of angular momentum L. The angular momentum L/N and the depletion of 
the ground states l^^i) are almost equal, depending on the fluctuations of the state. The dotted 
unconnected lines are the critical s-depletions, as estimated from MF theory; the calculations 
here are done for the permanents \nx,n,n, n), built over four orbitals, with equal occupations of 
the p-orbitals. The second lowest continuous line (yellow) determines the expectation value of 
the angular momentum Lmf/N, over (generally fragmented) MF states, which is given by the 
depletion |(1 — pi/N), here p\ = n\. The lowest continuous line (red) on the diagram depicts the 
maximum number of bosons N~* that can be loaded in a Gross-Pitaevskii condensate (marked 
as GP in the graph) without collapsing. For this curve one has to identify the axis 1 — ^ with 

N GP 

1 ^7— , i.e., p\ plays the role of the critical GP particle number at a given |Ao|(iV — 1). The 

difference in the estimation of the factor A cr = \Xo tCr \(N — 1) where the L = ground state collapses 

is evident; the overestimated value of A cr from the GP approach is larger than the MB one and 

puts the latter closer to the experimentally measured one Q]. The depicted stability behaviour 

does not significantly depend on the particle number N, making this behaviour universal (see text 

for more details). All quantities are dimensionless. 
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